

##########################################################################
# Title: Democratization and Religious Nationalist Mobilization Against ##
# a Small Minority: evidence from Myanmar                               ##
# Author: Megan Ryan                                                    ##
# Content: Maps and Main Analysis                                       ##
# Date: September24, 2025                                               ##                                                                ##
##########################################################################

#Clean the working environment and set the working directory

rm(list = ls())
setwd("C:/Users/memar/University of Michigan Dropbox/Megan Ryan/JMP") #Change to your working directory

##load required packages
library(sf)
library(dplyr)
library(ggplot2)
library(stargazer)
library(plyr)
library(lubridate)
library(readxl)


#read in datasets
Mil.ceremonies <- read_excel("Mil.ceremonies.xlsx")
Pcodes <- read_excel('Pcodes.xlsx') 
MBT <- read_excel('MBT.xlsx')  
MMR <- read_sf(dsn = "C:/Users/memar/University of Michigan Dropbox/Megan Ryan/JMP", layer = "mmr_polbnda_adm3_250k_mimu_1") #Change to your working directory

##Merge datasets

#select necessary columns
Pcodes <- Pcodes %>% select(SR_Pcode, SR_Name_Eng, Tsp_Pcode,Township_Name_Eng)

#Merge MBT dataset with township codes
#Mil.ceremonies_pcode <- left_join(Mil.ceremonies, Pcodes, by = c("Township_Name_Eng"), relationship = "many-to-many")
Mil.ceremonies_pcode.full <- full_join(Mil.ceremonies, Pcodes, by = c("Township_Name_Eng"))


#drop mismatched townships
Mil.ceremonies_pcode.full <- Mil.ceremonies_pcode.full %>% filter(!row_number() %in% c(260, 390, 408))

#sum total military ceremonies by township

colnames(Mil.ceremonies_pcode.full)[7] <- "TS_PCODE"
Mil.ceremonies_pcode.total <- Mil.ceremonies_pcode.full %>% 
  mutate(Total = if_else(is.na(SR_Name_Eng.x) == TRUE, 0, 1))
Mil.ceremonies_agg <- ddply(Mil.ceremonies_pcode.total,"TS_PCODE", numcolwise(sum)) 


#merge MBT activities dataset with Pcodes

MBT_pcode <- full_join(MBT, Pcodes, by = "Township_Name_Eng")


#drop mismatched townships
MBT_pcode <- MBT_pcode %>% filter(!row_number() %in% c(19, 176))


#sum total MBT activities by township
colnames(MBT_pcode)[8] <- "TS_PCODE"

MBT_pcode.total <- MBT_pcode %>% 
  mutate(Total = if_else(is.na(Region.state) == TRUE, 0, 1))
mbt_agg <- ddply(MBT_pcode.total,"TS_PCODE", numcolwise(sum)) 




##make plots
#create Myanmar shape files for plots
Mil.ceremonies_agg_geo <- left_join(MMR, Mil.ceremonies_agg, by = "TS_PCODE")
mbt.agg.geo <- left_join(MMR, mbt_agg, by = "TS_PCODE")

#only keep necessary variables

colnames(Mil.ceremonies_agg_geo)
Mil.ceremonies_agg_geo <- Mil.ceremonies_agg_geo %>% select(ST, ST_PCODE, DT, DT_PCODE, TS, TS_PCODE, TS_MMR, geometry, Total)
#only include desired variables
mbt_agg <- mbt_agg %>% select(TS_PCODE, Total)


#Figure 3: Map of MaBaTha Public Events by Township, 2013 - 2015

ggplot(mbt.agg.geo) + 
  geom_sf(aes(fill = Total))+
  scale_fill_viridis_c(option = "turbo")

#Figure 4: Map of State-Sponsored Buddhist Ceremonies by Township, 2006-2011

ggplot(Mil.ceremonies_agg_geo) + 
  geom_sf(aes(fill = Total))+
  scale_fill_viridis_c(option = "turbo")


#Merge both MBT and Mil.donations datasets for main analysis
Mil.ceremonies_MBT <- left_join(Mil.ceremonies_agg, mbt_agg, by = "TS_PCODE")



#add pro-military support variable (dichotomous variable whether the pro-military party (NUP)
#won the 1990 elections at the township level)

Mil.ceremonies_MBT <-Mil.ceremonies_MBT %>% mutate (NUP = ifelse(TS_PCODE == "MMR017017"|
                                                                         TS_PCODE == "MMR017019"|
                                                                         TS_PCODE == "MMR017001"|
                                                                         TS_PCODE == "MMR005037"|
                                                                         TS_PCODE == "MMR015021"|
                                                                         TS_PCODE == "MMR004008"|
                                                                         TS_PCODE == "MMR002006"|
                                                                         TS_PCODE == "MMR002005"|
                                                                         TS_PCODE == "MMR001017"|
                                                                         TS_PCODE == "MMR001006", 1, 0))


#add district variable for fixed effects

Mil.ceremonies_MBT.dt <- left_join(Mil.ceremonies_MBT, MMR, by = "TS_PCODE")

#add districts for townships not included in shape file

Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015201", "MMR015D221", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015202", "MMR015D221", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015203", "MMR015D221", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015301", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015302", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015303", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015304", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015305", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015306", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt<- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015307", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015308", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015309", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015310", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015311", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015312", "MMR015D331", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015313", "MMR015D332", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015314", "MMR015D332", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015315", "MMR015D332", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015316", "MMR015D332", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015317", "MMR015D332", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR015318", "MMR015D332", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR016319", "MMR016D333", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR016320", "MMR016D333", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR016321", "MMR016D333", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR016322", "MMR016D333", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR016323", "MMR016D333", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR013008", "MMR013D001", DT_PCODE))
Mil.ceremonies_MBT.dt <- Mil.ceremonies_MBT.dt %>% 
  mutate(DT_PCODE = if_else(TS_PCODE == "MMR013045", "MMR013D001", DT_PCODE))


## GLM Poisson models

poisson.model1 <- glm(Total.x ~ Total.y, family = poisson(link = "log"), data = Mil.ceremonies_MBT.dt)
poisson.model2 <- glm(Total.x ~ Total.y + NUP, family = poisson(link = "log"), data = Mil.ceremonies_MBT.dt)
poisson.model3 <- glm(Total.x ~ Total.y + DT_PCODE + NUP, family = poisson(link = "log"), data = Mil.ceremonies_MBT.dt)

summary(poisson.model1)
summary(poisson.model2)
summary(poisson.model3)


##Present results in table form

stargazer(poisson.model1, poisson.model2,poisson.model3, omit = 'DT', add.lines=list(c('Fixed effects','No', 'No', 'Yes')))


